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O , Abstract. The N = 3 Dicke model couples three qubits to a single radiation 

^^ , mode via dipole interaction and constitutes the simplest quantum-optical system 

^^ , allowing for Grcenberger-Horne-Zeilinger states. In contrast to the case N = 1 

Qv ■ (the Rabi model), it is non-integrable if the counter-rotating terms are included. 

The spectrum is determined analytically, employing the singularity structure of 
an associated differential equation. While quasi-exact eigenstates known from the 
C^ | Rabi model do not exist, a novel type of spectral degeneracy becomes possible 

Qj, which is not associated with a symmetry of the system. 
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1. Introduction 

The simplest model to describe light-matter interaction is the quantum Rabi model, in 
which a two-level system (two states of a single atom in the early applications) interacts 
with a single mode of the radiation field [f ] . A seminal step in its analytical treatment 
has been taken by Jaynes and Cummings through the invention of the "rotating- wave 
approximation" (RWA), valid close to resonance and for coupling strengths typical 
for atom optics [2]. The ensuing model can be solved analytically in a very simple 
way, because the RWA introduces a strong continuous symmetry [3], rendering it 
superintegrable [4] . A natural generalization of the Rabi model is the Dicke model, in 
which the radiation mode couples simultaneously to N two-level systems (qubits) [5] . 
It was first studied in the limit of large N, because it exhibits a phase transition to a 
"super-radiant" state for strong coupling [6, 7, 8]. Although the transition cannot be 
observed within atom optics [10], it should be realizable within circuit QED [11, 12, 13] 
and its equivalent has been experimentally observed in a Bose-Einstein condensate 
coupled to an optical cavity [14]. 

While these developments concern the Dicke model for large N, applications 
to quantum information technology have renewed the interest in the case of small 
N [15, 16, 17, 18, 19]. The model with three qubits allows in principle the 
dynamical generation of Greenberger-Horne-Zeilinger (GHZ) states [20] which could 
be of importance for future applications e.g. in quantum cryptography [19]. Possible 
realizations of the N = 3 Dicke model within circuit QED will be able to explore the 
strong coupling region [21] where the RWA is not feasible and one has to consider the 
full model. The [/(l)-symmetry induced by the RWA is so powerful that the Dicke 
modelj becomes integrable for arbitrary N [22], while the model including counter- 
rotating terms is non-integrablc for all N > 2 according to the criterion introduced in 
[4]. The case N = 1 is the only one where Schweber's technique [23, 24] or operator 
methods [25, 26] are applicable, therefore we shall employ in the following the method 
based on analysis of the associated differential equation in the complex domain [4] . 

The Dicke model for N = 3 is described by the Hamiltonian (h = 1), 

H' D = U aU+f£af + ^=(a + a^J2of. (1) 

i— 1 v i— 1 

Here u> denotes the frequency of the radiation mode, a (a^) is its annihilation (creation) 
operator, ujq is the energy splitting of the three qubits, described by Pauli matrices 
o~l' x 1 which are coupled through a dipole term with strength g' to the field. Because 
the qubits are equivalent, Hjj is rotationally invariant, leading to a splitting of the 
eight-dimensional spin-space into irreducible components, according to 

The N — 3 Dicke model is equivalent to two Rabi models and a system with spin 
S = 3/2. We shall confine ourselves in the following to the S = 3/2 model with 
four-dimensional spin-space. The Hamiltonian reads with w = 1, A = u>o/2 and 

g = g'/Vs, 

H D = a^a + 2AJ z + 2g(a + a i <)J x , (3) 

and J z and J x are generators of SU(2) in the spin-3/2 representation. Hjj possesses a 
Z2-symmetry (parity), P — e t1Ta a R, with the involution R acting in spin space as 

J Properly called "Tavis-Cummings model" if the RWA is used. 
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RJ Z R = J z , RJ X R --J x . We have PH D P = H D . The total Hilbert space L 2 (R)<g>C 4 
splits into two invariant subspaces (parity chains) labeled by the eigenvalues ±1 of P 
[21]. This discrete symmetry is familiar from the Rabi model and renders it integrable 
because P has as many irreducible representations as the dimension of the state space 
of the (single) qubit [4]. The same symmetry is present in the S = 3/2 Dicke model. 
However, it does not lead to integrability because a spin-3/2 has a four-dimensional 
state space and the eigenstates cannot be labeled by a quantum number representing 
the continuous degree of freedom (the radiation field) together with a second label for 
the discrete degree of freedom which would have to take four different values. The 
Z2-symmetry, providing only two different labels, is therefore too weak to make the 
S = 3/2 Dicke model integrable. Nevertheless, the symmetry leads to a considerable 
simplification of the analytical solution. 

2. The spectrum 

The operator R becomes a simple reflection after transformation of Ho into "spin- 
boson" form, using a SO (4)-transformation OHdO = H s b, with 

H sb = a i a + 2AJ x + 2g(a + a r )J z . (4) 

As in the Rabi model, the parity invariance can be used to partially diagonalize H s b. 
Define f ■■ ' ' ' ' 



1 
u =-f= 



/ 1 1 \ 

10 1 



(5) 



T -T 

V f o -f o / 

Then U'H s bU = H + + H_, where H± acts in H±; H+ and %_ are the two mutually 
orthogonal subspaces with fixed parity. We have 

H —'"^{% ± fr)-'' (o + <,,) (o ?)■ <6) 

We shall now represent the continuous degree of freedom in the Bargmann space £>, 
spanned by analytic functions f(z) [27]; T-L± is isomorphic to B<8)<C 2 . The operator T 
acts on elements of B as (Tf)(z) = f(—z). The eigenvalue equation H + ip = Eip takes 
with (p = (4>i(z), (f>2(z)) T the form of a non-local system of linear ordinary differential 
equations in the complex domain, 

z-^-Mz) + V5A<h(x) - Sgz^iz) - 3g— fa(z) = Efa(z), (7) 

az dz 

z-rfo{z) + >/3A0i(«) - gz<h(z) - g-rfciz) + 2A0 2 (-z) = E<j> 2 (z). (8) 

dz dz 

With the definitions (j>j(z) = cpj(—z), j = 1,2 and denoting the derivative with a 
prime, we obtain the following local system of the first order, 

(z-3g)<f/ 1 = (E + 3gz)<f>i-\/3A4>2, (9) 

(z-g)cf>' 2 = (E + gz)(f) 2 - VZAfa - 2Afo 7 (10) 

{z + 3g)& = (E-3gz)fa-T/3Afa, (11) 

(z + g)4>' 2 =(E-gz)4> 2 -V3A4> 1 -2A4> 2 . (12) 
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Figure 1. The singularity structure of system (9)— (12) in the complex plane. 
Solutions in powers of (z — z,) converge in the disks Dj with radius 1g around 
the regular singular points Zj = jg. The two ordinary points used to define the 
analytic continuations are z' = and zq = 2g. 



The system (9)-(12) has four regular singular points at z — ±<?, ±3g, and an irregular 
singular point at infinity [28] . The latter has rank 1 as in the Rabi model [29] , therefore 
a solution ip of (7,8) is an eigenvector of H + with eigenvalue E if and only if (fii(z) 
and 4>2{z) are analytic in the whole complex plane [27]. The indicial analysis of (9)- 
(12) shows that the exponents at the points z — ±3g are and E + 9<? 2 , whereas 
at z = ±g they are and E + g 2 . The exponent zero is three- fold degenerate at all 
regular singularities. This is the major difference to the Rabi model, which has only 
two regular singular points at ±g and the exponent zero is non-degenerate at each 
of these points. Formal solutions of (9)-(12) in terms of power series in z — Zj are 
possible in regions Dj with Zj = jg, j = ±1,±3, see Fig. 1. The expansion around 
Z\ = g reads for j = 1, 2, 



4>j{z) 



oo 

£< 

71=0 



-,n(z-g) n , 0j(z) 



J2*jAz~9T- (13) 

The expansion (13) is absolutely convergent in _D 1; with radius of convergence 2g and 
4>j{z), 4>j{z) are analytic at Z\. Likewise, there is an expansion of <ftj(z), <f)j(z) around 
z 3 , which we denote as ipj(z),ij>j(z), 



ipj(z) 



n=0 



,(z-3g) n , $Az) 



n=0 



3g)". 



(14) 



(14) converges in Dj, and ipji z )i 4>j( z ) are analytic at z%. Using the identification 
4>j(z) = <f>j{—z) and ipj{z) — ij}j(—z), (13) and (14) furnish series expansions of <pj( z ) 
in regions D-\ and -D-3, respectively. These series lead for arbitrary E to functions 
which are analytic at their expansion points but develop branch-cuts at the other 
singular points. The discrete set of eigenvalues E n , n = 0, 1,2, ... is determined by 
the condition that all four expansions describe the same function </>_,- (z), i.e. that they 
are analytic continuations of each other [29] . 

Inserting the ansatz (13) into (9)-(12) yields the following coupled recurrence 
relations for n > and with x = E + g 2 , 

2g(n +l)(x- (n + l))a 2 ,„+i = (3A 2 - 



+ 2xn - 2g 2 



2xg 2 )a 2 . 
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+g{\ -Ax- 2g 2 + An - 4)a 2 ,„_i - 3.g 2 a 2j „_ 2 (15) 

+AAg(n + l)a 2 ,„+i + 2A(x + 2g 2 - n)5 2 ,„ + 6Aga 2 , n -i, 

and 

8 5 2 (n + 2)(n + l)a 2: „ +2 = -2g(n +l)(3n + 2-3x + 8g 2 )a 2 . n+1 

+ (3A 2 - n 2 + n(2x - 16.g 2 ) - Ag 2 - (x - 2g 2 )(x - Ag 2 ))a 2 . n 

+ (Ax - 10g 2 + 3 - An)a 2 , n -i - 3g 2 a 2 ^ 2 (16) 

-8Ag(n + l)a 2 ,„+i + 2A(x - 4g 2 - n)a 2 ^ n - QAga 2 , n -i- 

The coefficients ai, n , ai in read in terms of the a 2 . n , a 2 . n as 

ai, n = -=— [(n - a;)a 2jn -.ga 2 , n -i + 2Aa 2j „], (17) 

ai,„ = 7f7"[( n - » + 2g 2 )a 2 ,„ + ffa 2l „_i + 2g(n + l)a 2 ,„+i + 2Aa 2) „]. (18) 

The system (15,16) cannot be reduced to a linear three-term recurrence relation, 
therefore continued- fraction techniques are not applicable to the present model [30]. 
Moreover, the relations (15)— (18) are not sufficient to determine the eigenfunctions 
uniquely because there are three linear independent solutions of (15,16) which are 
analytic at z\. The initial conditions for (15,16) are given by the set {a 2 ,o, a 2 ,o, ^2,1} 
and the three solutions {^ '(z),(j>j (z)} of (9)-(12) analytic at z — g are obtained 
by setting one element of the set to 1 and the others to 0. 
The general solution of this type is therefore given as 

3 3 

&(«;{7fc» = I>fc^ fc) (*), &(*;{7fc}) = X>fc$ fc) (*0, (19) 

fe=l k=l 

with 7fe to be determined. In the same way we obtain for the ^pj{z) 1 ^j{z) the 
recurrences, 

2g(n + l)(n-x- 8.g 2 )a 2 ,„+i = (3A 2 - n 2 + n(2x + 16.g 2 ) - (x + 8g 2 )(x + 2g 2 ))a 2 . n 
+g(An -3- Ax- 14g 2 )a 2 ,„_i - 3.g 2 a 2i „_ 2 (20) 

+2A(a; + 8g 2 — n)a 2 ^ n + 6Aga 2 ,n-i, 

and 

24.g 2 (n + 2)(n + l)a 2j „ +2 = 2g(n + l)(5x - hn - 3 - 32g 2 )a 2 ,„ + i 

+ (3A 2 - n 2 + n(2x - 32g 2 ) - 6.g 2 - (a; - 10.g 2 )(a; - 4 5 2 ))a 2 ,„ 
+g(Ax - An + 3 - 22.g 2 )a 2 ,„_i - 3.g 2 a 2: „_ 2 (21) 

-12A.g(n + l)a 2 ,„+i + 2A(x - IO3 2 - n)a 2: „ - 6A#a 2 ,„_i. 

For the <l\, n , a\. n we have, 

ai,„ = -=— [(n - x - 2g 2 )a 2 . n - ga,2, n -i + 2g(n + l)a 2 , n+ i + 2Aa 2 , n ], (22) 

ai,„ = — 7f7"[( n - x + 4 5 2 )«2,™ + ga,2, n -i + 4g{n + l)a 2 , n+ i + 2Aa 2 , n ]. (23) 

The solution space of (9)-(12) analytic at z = 3g is likewise three-dimensional, 
determined by the initial set {a 2i o, a 2 . 07^2,1} and the general solution reads 

3 3 

^(z;{c k }) = Y,Ck^ k) (z), i> j (z;{c k }) = Y / Ck$ k) (z), (24) 

fe=i fe=i 
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with three unknown constants c^. (19) and (24) are the most general solutions analytic 
at g and 3g respectively, if the spectral parameter x is not a positive integer (see (15)) 
or satisfies x + 8g 2 = n, n = 0, 1,2, . . . (20). These special values for x determine the 
two types of baselines in the model where the exceptional spectrum is located (the 
Rabi model has only one type of baseline). We shall now determine the condition 
under which all four sets of series expansions for <j>i(z),<j)2(z) describe the same 
functions, which are therefore analytic in the whole complex plane and correspond 
to an eigenvector (<j)i,<j)2) T of H+. 

Because the functions {<f>i(z),(f>2(z), 4>i(z), 4>2(z)} satisfy the same differential 
equation of the first order as {ipi(z), tp2(z), 'ipi(z) , ^(z)} , both sets will coincide in 
D2 = Di n Ds, if they coincide at one regular point zq E D2 [29] . This yields four 
equations for the functions in (19) and (24). Furthermore, {<pi{z),<j>2{z),^i(z),^2(z)} 
and {0i (— z),4>2{— z),<f>i(— z),4>2(— z)} satisfy the same differential equation and 
coincide in all of Do = D\ n D_\ if they do so at a point z' e Do. Obviously, only 
two of the four equations are independent if z' — 0. (4>i(— z), (j>2{— z)) T is then the 
analytic continuation of (</>i(z), 4>2(z)) T into the disk D-\. But because 4>j(z) = i>j(z) 
for z € D2, it follows that ipj(-z) is the analytic continuation of 4>j(—z) (and therefore 
of 4>j{z)) into the disk D-3. The six equations 

(j)j(zo) = i>j(z ), <j>j(z ) = i>j(z ), (f>j(0) = 4>j(0), (25) 

for j = 1, 2 and Zq € D 2 are equivalent to the analyticity of 4>j{z) in C. A non-trivial 
solution of (25) can be found if the parameters {7^, c^} are not all zero. The functions 
4>j(z), <f>j(z), . . . depend parametrically on the energy E = x — g 2 . Define then the 
matrix M + {x, zq) as 

( ^\z ) 4 2) (z ) 4 3) (^o) 

^\z ) ^\z ) 4 3 \z ) 

^\z ) ^\z ) ^\z ) 

$\zo) J2 2 \zo) ^ 3) (zo) 



V 

with <5<^ fc) (0) = 4>f } (0) - 4>f } (0). It follows that the G-function of the Dickc model 
for positive parity, 

G+(x,zo) = detM + (x,z ), (27) 

is zero for all zq G D2, if and only if a; = E + g 2 corresponds to an element 
of the spectrum of H + . The discrete set of zeros x n with G+(x„,zo) = for 
n = 0,1,2,... determines the regular spectrum a r (H + ) — {x n — g 2 } n eTN a of H+ 
[4]. The regular spectrum of _ff_ is given in an analogous manner, starting from 
recurrences (15,16,20,21) by reversing the sign of the coupling terms between a n and 
a n , resp. a n and a n (see (6)), constructing functions (jr- _(x, z), m _(x, z), . . . and the 
matrix M_(x, zq). Fig. 2 shows G±(x, 2g) for g = 0.25 and A = 0.7. 

The functions m _ (x, z), (fn _ (x, z) entering M_ (x, zq) are not related in a simple 
manner to their counterparts (f)j + (x, z), 4>^ + (x, z) for positive parity, in contrast to the 
case of the Rabi model. There one has only one pair of functions {<fi±(x, z), 4>±(x, z)} 
and the G-function of the Rabi model reads simply 

G%(x) = $±{x,0) - <f> ± (x,0). (28) 



M + = 



-tfHzo) 


-4> { i\zo) 


-<t> ( i\zo) \ 


-<&\zo) 


-4> ( 2 2 \z ) 


-#(^0) 


~^\zo) 


~4>f\zo) 


-^\z ) 


-^{ZO) 


-$%) 


-$\z ) 


5$\o) 


^i 2) (0) 


^i 3) (0) 


6$\o) 


S^iO) 


^ 3) (o) / 



(26) 
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Figure 2. The G- functions of the Dicke model for positive (red) and negative 
(blue) parity. A = 0.7 and g = 0.25, therefore the baselines of first kind are 
located at integer, and the baselines of second kind at half-integer values. At 
both kinds has G±(x,0.5) poles of order 1 or 3. 



Moreover, 4>-{x,z) — <f>+(x,z), <f}—(x,z) = —4>+{x,z). It follows at once that the 
regular spectrum of the Rabi model is not degenerate between states of different parity, 
as G+(x) = G^(x) = implies c/>±(x,z) = 0. On the other hand, regular states with 
the same parity cannot be degenerate because the formal solution <f>+(x,z), analytic 
at Z = .9, is unique for x ^ IN. The only possibility for degenerate eigenvalues in 
the Rabi model occurs therefore in the exceptional spectrum, where x £ IN. The 
two degenerate states at these points have different parity because there are only two 
linear independent formal solutions of the eigenvalue equation and the pole of G± (x) 
at integer x is lifted in both G5(x) and G:?(x), which corresponds to a condition 
satisfied by the model parameters g and A [31, 4]. 

As mentioned above, G±(x, zq) has two different types of baselines, located at 
x = n, n £ IN (first kind) and x — n — 8g 2 , n £ IN (second kind). Although G±(x, z ) 
has pole singularities in x at these values for general g, A, there may be (exceptional) 
eigenvalues E\ — n~ g 2 , resp. E 2 — n — 9g 2 , if the singularity is lifted in G+(x, z ) 
or G^(x,zo). This exceptional solution, however, is usually non-degenerate, because 
there is no single lifting condition (as in the Rabi model), valid for both parities. § 
Therefore, a quasi-exact spectrum in the sense of the Rabi model is not present in the 
Dicke model. If one defines the quasi-exact spectrum differently, by demanding that 
the eigenfunctions are polynomial in z (apart from a common factor), this possibility 
is not ruled out in principle, although the set {n,g, A} would then have to satisfy 
three consistency equations [31], making a solution with integer n unlikely. In fact, as 
these solutions lie necessarily on baselines and are parity degenerate, the consistency 
equations given by Kus and Lewenstein must comprise the two independent lifting 
conditions for G±(x, Zq). On the other hand, G+(x, z ) — G®(x,Zo) is no longer 
tantamount to vanishing of the wave-functions (j>j,±(x, z) themselves, therefore regular 

§ Non-degenerate exceptional solutions appear in the Rabi model as well, but then the pole is not 
lifted in <f>±(x,0) but only in the difference, either 0_|_(a;,O) — </>+(a;,0) or rp—(x,0) — <f>—(x,0). These 
solutions are exceptional but not quasi-exact. 
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Figure 3. The spectrum of the Dicke model for even (red) and odd (blue) parity 
at A = 0.7 and for varying g. The j/-axis shows x = E + g 2 , baselines of first 
kind arc horizontal straight lines. The ground state has odd parity as in the Rabi 
model. The two ladders of eigenvalues with different parity intersect within the 
regular spectrum. There are no degeneracies (but narrow avoided crossings) for 
fixed parity in this parameter window. The baselines of first kind (not depicted) 
and of second kind (dashed lines) emerge as limiting values in the deep strong 
coupling regime g > 1. 




Figure 4. The Rabi spectrum for the same parameters as in Fig. 3. Degeneracies 
occur solely between states of different parity and are always located on the 
baselines. 



eigenvalues of the Dicke model may well be parity degenerate. Figs. 3 and 4 show the 
Dicke and Rabi spectra for fixed A and varying g. It is apparent that the coupling 
between exceptional eigenvalues and degeneracies renders the Rabi spectrum much 
more "regular" than the Dicke spectrum, apart from the complication coming from 
two kinds of baselines in the latter. 

All regular eigenvalues of the Rabi model with fixed parity correspond to unique 
cigenfunctions because the exponent of the indicial equation at each singular point 
is non-degenerate [28, 31] and there exists at most one solution analytic at z = ±g. 
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In contrast, we have three solutions analytic at each of the Zj in the case of the 
N = 3 model. Although generally there is only one solution analytic at all four 
singular points, the possibility is not excluded that the kernel of M±(x, z ; <?, A) has 
dimension > 1 at some value x, which would correspond to a degeneracy within 
a given parity chain H±. Indeed, for this it is necessary that the linear term in 
the characteristic polynomial of M±(x, zo] <?, A) vanishes, providing a condition to be 
satisfied by the model parameters g and A, in analogy to the equation determining 
the quasi-exact spectrum of the Rabi model [31]. On has to note here that this 
equation is not independent from det M±(x, Zo;g, A) = 0, because the value of x is 
not restricted to integers. Both equations form a coupled system to determine the 
triple (xdcg,3dcg, Adcg) where two states with equal parity are degenerate. If such a 
triple exists, a level crossing at E — x^ cg — g^ will appear within the corresponding 
parity chain H±. These degeneracies would not originate from a global symmetry of 
the model and thus are not accidental, because the degenerate states do not belong 
to dynamically decoupled subspaces. The numerical investigations done so far have 
shown no hint to this novel type of degeneracy yet. However, in a recent work on 
the N — 2 model with inequivalent qubits, Chilingaryan and Rodriguez-Lara have 
discovered level crossings within spectra with fixed parity [32]. This is a strong 
indication that the phenomenon predicted here for the N = 3 model is not forbidden 
by some special feature of the matrix M±(x, zq; g, A). 

3. Conclusions 

We have computed the spectrum of the Dicke model for three qubits analytically using 
the technique based on formal solutions in the Bargmann space, where the spectral 
condition corresponds to analyticity in the whole complex plane [27, 23, 4]. In contrast 
to the N = 1 model, the argument used by Schweber to derive a continued-fraction 
representation of the spectral condition is not applicable, because the formal solutions 
are not given in terms of linear three-term recurrence relations [30]. If one defines 
the model on a truncated Hilbcrt space, matrix-valued continued fractions could be 
employed in principle [33], but this approach suffers from ambiguities and can be 
justified only in the case N = 1, by its formal equivalence to Schwebcrs method [34]. 
The Dicke model shows much more spectral "irregularities" than the Rabi model, 
whose spectral graph is restricted by the position of the quasi-exact eigenvalues, 
confining degeneracies to the baselines. The N = 3 model possesses two types of 
baselines, but they govern only the asymptotics for strong coupling (Fig. 3); the quasi- 
exact spectrum does not exist and levels corresponding to different parity intersect 
within the regular spectrum, while the exceptional spectrum is non-degenerate. On 
the other hand, the structure of the G-function for N = 3 predicts degeneracies 
within the parity chains, which are forbidden for N — 1. These degeneracies are not 
related to a symmetry of the model and do not imply intcgrability. They cannot be 
termed "accidental" either, as the degenerate states do not belong to different invariant 
subspaces. However, if they appear in sufficient number, which seems to be possible 
for large N, it could lead to a level statistics resembling Poissonian behavior, which 
has been found numerically for N > 20 in the coupling region below the quantum 
phase transition [35] . The implications of this novel type of degeneracy for the notion 
of (non)-integrability in systems with less than two continuous degrees of freedom [4] 
have yet to be explored. 
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